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ABSTRACT 

The dark matter halos in ACDM cosmological simulations are triaxial and highly flattened. In many 
cases, these triaxial equilibria are also tumbling slowly, typically about their short axes, with periods 
of order a Hubble time. Halos may therefore exert a slowly changing external torque on spiral galaxies 
that can affect their dynamical evolution in interesting ways. We examine the effect of the external 
torques exerted by a tumbling quadrupolar tidal field on the evolution of spiral galaxies using N-body 
simulations with realistic, disk galaxy models. We measure the amplitude of the external quadrupole 
moments of dark halos in cosmological simulations and use these to force disk galaxy models in a series 
of N-body experiments for a range of pattern speeds. We find that the torques are strong enough 
to induce long lived transient warps in disks similar to those observed in real spirals and also induce 
the bar instability at later times in some galaxy models that are otherwise stable for long periods of 
time in isolation. We also observe forced spiral structure near the edge of the disk where normally 
self-gravity is too weak to be responsible for such structure. This overlooked influence of dark halos 
may well be responsible for many of the peculiar aspects of disk galaxy dynamics. 

Subject headings: methods: n-body simulations - methods: numerical -cosmology: dark matter halo 



1. INTRODUCTION 

The nearly flat rotation curves of spiral galaxies 
are direct evidence for the existence of massive dark 
halos around the galaxies, wi thin the paradigm of 
Newtonian-Einsteinian gravity (IRoberts fc Whitehurstl 
19751: iRubin et all fl979l [1983 119851: Ivan Albada et all 
198a iSalucci fc Burkertl 12000: ISofue fc Rubinl l200if l 
The dark halos forming in cold dark matter cosmolog- 
ical models have been identified with the inferred dark 
halos of galaxies; a substantial amount of work has gone 
into trying to tes t for the consistency between reality and 
simulation (e.g., iBarnes fc Efstathioul 119871: iFrenk et al.1 
19881: iDubinski fc Carlberd 119911: I Navarro et all 119961: 
Bullock et all 120011: iNeto et alj I2007D . Most work has 
focused on the spherically-averaged density profile and 
the implications of this profile for galactic rotation 
curves. Although debate still con tinues on the value 
of the slope of th e inner cusp (jGentile et al.l l2004t 
iFerreras et al.|[2b07f ) and whether or not th e inferred ro- 
tation curves are consistent with real ones (lMoorolll994l 
Ide Blok fc McGaughlll998l : Ide Blok et al.ll2003h . it seems 
that the consistency between the real and simulated 
dark halos is quite strong and that the cold dark mat- 
ter (CDM) par adigm is in reasonable shape in this sense 
(|Primackll2007l ). 

Other work has focused on secondary dynamical 
effects of dark halos that arise from their triaxial 
nature and their implication for other properties in 
disks, namely, oval distortion s , warps and perhaps bars 
(iBinnevI fl97& ISparkd 11984 iFranx fc de Zeeuwl 119921: 
Sackett et al.lll994l:lKuiiken fc Tremaindll994l : IWeinberg| 
19981: iSellwoodl 120031: iBerentzen et al.ll2006D . The flat- 
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tened potential of triaxial halos has been recognized in 
some work to induce oval distortions to otherwise cir- 
cular disks, as well a s vary the velocity along the orbit 
(jHayashi et alJ 120071 ). The orbits of stars in the disk 
become flattened in a direction orthogonal to the halo 
density. However, the detection of such oval distortions 
and non-circular velocities is marginal at best, suggest- 
ing that halos may be nearly axisymmetric in the plane 
of an ordinary, high surface-brightness spiral galaxy. On 
the other hand, the effects of non-circularity in dwarfs 
and low-surface brightness galaxies may be s tronger 
(|Valenzuela et alj|2007t lHavashi fc Navarroll2006D . 

Another source of complex dynamical effects is the 
probable misalignment between a disk and a triaxial dark 
halo. Disk-halo misalignment has been argued to be the 
origin for warps seen in most edge-on spiral galaxies. A 
misalignment implies that the disk will feel a torque from 
the halo and the nearly circular orbits will therefore pro- 
cess. If the disk has self-gravity and the halo potential 
is static, warped m odes can arise (e.g., iToomrd 119831 : 
iSparke fc Casertanol IT9881 : lKuiikenlfl99lh . A fatal flaw 
in this hypothesis is that gravitational interactions be- 
tween the disk and the halo are strong, with the process- 
ing d isk experiencing dynamical friction from th e halo 
(e.g.. lNelson fc Tremainell 19951: iBinnev et al.lll998T ). One 
therefore expects the disk and halo to become aligned 
with one another within a few dyn amical times and this 
is ind eed borne out in experiments ([Dubinski fc Kuiikenl 
1995). Thus it has been suggested that it is the outer halo 
that is misaligned with the disk, while a ti ght coupling 
exists between the inner halo and the disk (jBailin et alJ 
127)051 IBinnevI 12001 . 

The origin of galactic warps then remains a perplexing 
issue but alternative scenarios suggest a transient origin. 
Other ideas suggest that the cosmic infall of gas and 
dark matter alters the relative orientation the disk and 
dark halo due to addition of angular momentum (e.g. , 
IQstriker fc Binnevl 119891 : iDebattista fc Sellwoodl 119991 : 



2 



I Jiang fc Binnevl 1 1999ft and th e resulting torques on the 
disk lead to transient warping. Ivan der Kruitl (|2007ft also 
considers the possibility of the onset of a late warped 
outer disk, due to gas infall, in a configuration that is in- 
dependent of the relatively younger inner truncated stel- 
lar dislc^^h^JnMl picture has also been recently used 
bv lShen fc Sellwoodl (|2006ft to reproduce the warping of 
a simulated disk galaxy; in fact, their simulations corrob- 
orate the result that the line-of-nodc of the warp forms 
a leading spiral. 

Tidal fields from satellite galaxies are known t o excite 
dramatic behavior in spirals; iKalirai et alJ (|2006ft invoke 
the possibility of a warp or overdensity induced by satel- 
lite interaction to explain the origin of the secondary cold 
population that they identify in M31 while the grand de- 
sign spiral in M51 is inferred to be th e result of it's inter- 
action with its companion NGC 5194 (|Toomre fc Toomrd 
Il972t ISalo fc Laurikainenll2000ft The Magellanic origin 
of the Galactic warp (Weinberg 1998; Garcia-Ruiz et al. 
12002ft is revisited by iWeinberg fc Blitz! (|2006l ) in which 
they suggest that the warp is formed as a joint handi- 
work of the tidal field of the Magellanic Clouds as well as 
the effect of the distortions that such has on the Galactic 
halo. The overall picture that emerges from these studies 
is that the source of the warping torque lies with asym- 
metries in the halo, at radii beyond the gravitational 
influence of the disk. 

Dark halos formed in simulations performed within 
the CDM scenario suggest that not only do they often 
settle into triaxial figures of equilibriu m but that they 
can also tumble, much like a rigid body (iDubinskil Il992t 
iBureau et all Il999t iBailin fc Steinmetzl |20o1k Static 
and tumbling equilibrium are both perfectly valid in 
analogy to the two solutions for the homogeneous Ja- 
cobi ellipsoids (|Chandrasekharl 11969ft which include a 
tumbling and stati c solut ion with internal circulation. 
IBailin fc Steinmet3 (|2004ft find a log normal distribu- 
tion of halo pattern speeds with a mean value of about 
flp = 0.15 km s _1 kpc -1 , corresponding to a tumbling 
period of 44 Gyr. A typical halo will then tumble through 
120 degrees over a Hubble time. The tumbling periods 
of dark halos are rather long but nevertheless, a pic- 
ture emerges in which disk galaxies are embedded within 
slowly tumbling, triaxial halos with rotation axes that 
are probably slightly misaligned with the disk rotation 
axes. It is likely that the inner halo and disk are aligned 
with each other but the outer halo (say beyond 100 kpc) 
is tumbling slowly and thereby affecting the disk with its 
tidal field. A pertinent question to ask then is how strong 
are the torques that a tumbling dark halo exerts on the 
disk at the center and whether or not these torques can 
have a significant influence on the evolution of a spiral 
galaxy over a Hubble time? 

In this paper, we carry out a series of experiments to 
study the effect of a slowly tumbling, external quadrupo- 
lar potential, on the evolution of a disk, with the primary 
goal of inducing galactic warps. We first measure the ex- 
pected strength of the quadrupolar tidal field directly 
from dark matter halos within ACDM cosmological sim- 
ulations and subsequently estimate the torques expected 
on exponential disks. We find that the typical torques 
from cosmological dark halos are generally an order of 
magnitude larger than the effect of the LMC, rendering 
them more effective in trigerring something dynamically 



interesting. 

Guided by such experimentally motivated numbers, we 
set up simulations of ideal galactic models embedded 
in nearly spherical halos, forced slowly by an external 
quadrupolar tidal field. We discover in these experiments 
that it is fairly simple to make warps of amplitudes sim- 
ilar to the observed ones. 

We also find that the gradual tidal forcing of tumbling 
dark halos has the effect of pushing an otherwise stable 
disk into a region where it becomes subject to the bar 
instability at late times. This is a new bar triggering 
mechanism and may be an additional factor that affects 
the evolution of t he fraction of barred galaxies over c os- 
mic history (e.g.. ICurir et al.ll2008t TSheth et al.ll2008ft . 

The plan of the paper is as follows. In §2, we calculate 
the amplitude of the torque on a galactic disk, expected 
from the external tidal field of the triaxial dark halo, as 
extracted from cosmological A-body simulations. In §3, 
we present A-body models of M31 and the Milky- Way 
that are forced by the external quadrupolar field with the 
same strength expected from cosmological simulations. 
In §4, we present the results of these simulations while 
in §5, present analytical models of rigid disks forced by 
external tidal fields, to quantify these effects. In §6 we 
conclude with a discussion of the implications of these 
results for warps, bars, and spirals in real galaxies. 

2. TIDAL TORQUES ON DISKS FROM COSMOLOGICAL 
DARK HALOS 

The dark halos in CDM cosmological models are 
highly flattened and triaxial wit h typical axis ratios 
c/a = 0.4 and b/a = . 6 (e.g., iBarnes fc Efstathio ' 
1987t iFrenk et all Jl988t IDubinski fc Carlberd 119911 



Warren et all 11992 : Uing fc Sutdl2002l ). Analysis of the 



halo spin parameter shows a tendency for alignment with 
the minor axis, suggesting that a disk that forms within 
a dark halo is likely to have its own spin angular mo- 
mentum closely aligned with the halo minor axis (e.g., 
| Dubinskilll992t IWarren et al.lll992l : IDubinski fc Kuiikenl 
Il995t lBinnevll2007ft . The angular momentum within ha- 
los is distributed between internal streaming and tum- 
bling motion. Early studies showed that halos are slowly 
tumbling through space with long periods. More recent 
quantitative analyses in standard ACDM cosmological 
models show that the dis tribution of tumbling freq uen- 
cies is roughly log- normal (jBailin fc Steinmetzll2004ft and 
peaked at about 0.15 km s -1 kpc -1 with a = 0.83. The 
dissipative infall of the baryonic matter during galaxy 
formation can also modify t he halo, both increasing its 
central concentration (e.g., iBlumenthal et all 11984ft as 
well as rounding out the shape by increasing b/a within 
the region of influenc e of the forming disk while leaving, 
c/a a bout the same (|Dubinskil [1994; K azantzidis et all 
2004). The axis ratios of the volumetric density con- 
tours at radii beyond the influence of the disk are less 
affected and remain about the same as the initial values. 

Disks and halos are also found to be tightly coupled 
through dynamical friction in the inner regions of galax- 
ies (Dubinski & Kuijken 1995) so one would expect the 
disk to lie in the principal plane of the halo, proba- 
bly with the spin vectors of the disk a nd halo aligned 
(jBailin et al.ll2005HLibeskind et al.ll2007ft . However, mis- 
alignment may persist at larger radii and the slow tum- 
bling of a halo could lead to an external torque on 
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the disk. Although current self-consistent simulations 
of galaxy formation have not fully addressed this point, 
we make the hypothesis that the disks are aligned with 
the principal planes of their inner halos through dynam- 
ical friction while the outer halos remain misaligned and 
may slowly be changing orientation with the tumbling 
frequencies measured in cosmological dark matter simu- 
lations. 

We estimate this torque through the analysis of more 
than 2000 dark halos extracted from a cosmological 
dark matter simulation with near standa rd parameters 
ft A = 0.7, n m = 0.3, a 8 = 0.9 and h = 0.7 (|Spergel et all 
2007). The simulation contains 512 3 particles in a 
cube of dimension L = 70h~ Mpc and was run using 
the G OTPM cosmological N-body code (Dubinski et al. 
2004). The simulation was run from z = 70 using 
2800 equal time steps. The co-moving particle soften- 
ing length was set to 3.5/i _1 kpc. Halos are extracted in 
a two stage process; we use a friends-of-friends algorithm 
to determine the dense centers of candidate halos. These 
regions could either be associated with an isolated halo 
or merging pair or group. The particles in a spherical 
region surrounding these dense centers are extracted out 
to roughly the distance where the density drops to 200 
times the critical density and then analyzed separately. 
We fit Navarro, Frenk & White (1996, NFW) profiles 
to all these candidates but retain only those which have 
fit profiles that lie within a threshold x 2 - I n this way, 
the original sample of about 2700 halos is reduced to a 
sample of about 2200 halos that have well-defined pro- 
files, undisturbed by major mergers or significant asym- 
metries. We then take a careful look at this sample to 
learn something about halo shapes and torques. 

Many studies have determined th e distribu- 
tion of ellipsoidal shape o f hal os dFrenk et al 
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1991 [Cole fc Lacevl TlOM I.Iing fc Sutol l200l 
Bailin ~fc Steinmetzl l2005T h We take this analysis 



one step further and calculate the distribution in ex- 
pected strengths of tidal torques on a disk that may be 
misaligned with the outer halo. An explicit assumption 
is that the disk and halo are aligned and will remain 
so because of dynamical friction in the inner regions 
while misalignments between the outer halo will persist 
and change dynamically as the halo tumbles through 
space. We find that the relative orientations of the inner 
(r < r s ) and outer halo (r > r s ) is greater than 20° for 
half of the halos in our cosmological simulation (Fig. [1]). 
We therefore expect significant misalignment that can 
lead to torques on the disk from the outer halo. We 
proceed by analysing the potential of the outer halo 
using a standard multipole expansion method. 

We are interested in the possible torques at the cen- 
tre of a disk tilted with respect to the outer halo. Pre- 
vious work has shown that disks and inner halos are 
aligned within a radius roughly equal to the extent of 
the disk, though su c h ali g nment may extend further 
(|Dubinski fc Kuijkenl 119951 : iBinnev et all 1 1998t h Any 
torque acting on the disk from the halo probably then 
arises from misalignments beyond this radius. It is use- 
ful to quantify these effects more concretely. According 
to the NFW model, dark matter halos extend to the virial 
radius r2oo = cr s where c is the concentration, r s is the 
scale radius of the density profile and the mean enclosed 



density is 200 times the critical density p c = 3Hq/8ttG. 
At r2oo, halos have characteristic mass M200 and circu- 
lar velocity v% 00 = GM200/ '^200- Typical concentrations 
for g alactic scale halos are c w 15 (e.g., iBullock et al.1 
1 2 lh though they can vary considerably with a scatter 
A(logc) = 0.18. A galaxy halo is usually characterized 
by the peak value of the rotation curve v c . max rather 
than i?2oo- Analysis of the NFW rotation curve shows 
that the maximum value of the rotation curve v Ctmax oc- 
curs at r = 2.16r s and is related to W200 : 

W200 = 2.15 Vc.max /(c) (1) 

1/2 



with 



/(c)= C - 1 /2(l n (l + c )_ 



1 + c 



(2) 



Given the parameters c and v c<max we find i>2oo and then 
can deter mine r^no (or r„) and M200 through the usual 
identities (N avarro et al.lll996l ): 
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As a specific example, consider a Milky Way sized 
NFW dark halo with with a v Cimax — 220km s _1 and 
concentrations with the range c = 10 — 20. Assuming 
h = 0.7, these models will have virial masses in the 
range M200 = 1-1 — 2.0 x 10 12 M Q and scale radii from 
r s = 11 — 26 kpc. The scale radius r s is roughly the size 
of the disk and since the disk and halo are tightly coupled 
within this region we expect alignment. We can expect 
misalignment and tidal torques from the halo mass dis- 
tribution for halo mass beyond the edge of the disk and 
so we need to calculate the external component of the 
potential beyond r > r s to estimate the tidal torques. 

For a halo particle distribution, the external potential 
from matter with r > r at a point (r, 9, </>) with r < r 
is given by the expression (e.g., Binney and Tremaine 
1987), 



$ext(r;r < r ) = -4ttG 



l=oo m—L 

E E 

1=0 m=-l 



c lm r l Yr(0A), (5) 



where Y™ is the usual spherical harmonic function and 
the coefficients ci m are evaluated from the matter beyond 
r > 7*0 through: 



N -(7+1) 
Cim = > , Tl" , 7T *r*(fla, 4>a), T a > r (6) 



S mm + 1) 

with a indexing a list of the N particles that have r > r$. 

For a dark halo, it is natural to compute these coeffi- 
cients for the particles with r Q = r s and 2r s , as a reason- 
able measure of the external torquing potential of a dark 
halo. 

Taking the lead from CMB analysis, a useful way of 
quantifying the strength of the halo tidal potential is 
through the parameter q defined as: 



/ 171 — 1 
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\7n— — l 



1/2 



\Ch. 



21 + 1 
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The monopole term will lead to zero torque and the 
dipole terms are probably unimportant since dark halos 
tend to be ellipsoidal and do not show a significant lop- 
sided mode. The most important terms for torquing are 
the quadrupole terms 1 = 2; one expects successive even 
I terms to have some effect but with increasingly lower 
strengths. We therefore focus on the quadrupole terms 
that have a strength given by the parameter c 2 (r > ro) 
as the main source of torque on a misaligned disk. We 
now calculate this value for our sample of halos and de- 
termine its importance for disk torques. 

Since the value of c 2 depends on the mass of the halo, 
we can rescale all halos to some reference mass. We 
are interested mainly in the distribution of the relative 
strength of the external tidal potential in typical spiral 
galaxies; so we first renormalize the halos to have the 
same scale radius and scale mass. The following formal- 
ism is then undertaken with these halos. 

1. We first fit the scale parameters r s and M s to the 
NFW density profile in log p — log r space with the 
potential-density pair of the spherical NFW profile 
given by: 



$( r ) : 
p(r)-. 



GAL 



log(l + r/r s ) 



Att r(r + r s ) 2 



(8) 
(9) 



2. Next we determine r 20 o as the radius containing 
the mean density 200p c and find c = r 2 oo/r s and 
so compute v c<max: v 2 oo and M 200 = M s [ln(l + c) - 
c/ (1 + c)] from the potential. 

3. We rescale all halos to a putative Milky Way model 
for our determination and call the rescaled tidal 
parameter C2,mw ■ We use the parameters v c . max = 
220 km s -1 , a concentration of c = 10 which imply 
r s = 26 kpc and M 200 = 2 x 10 12 M Q from the 
equations above. 

To correlate this parameter to a specific example, we 
compute the expected value of the quadrupole tidal pa- 
rameter for a satellite of the type of the Large Magellanic 
Cloud, (mass M = 10 10 M at a distance of R = 50 
kpc). We refer to this tidal parameter as c 2 ,la/c- The 
LMC has a minor tidal effect on th e dynamics of the 
Gala x y at its current distance (e.g.,. iHunter fc Toomrel 
1969: G arcia-Ruiz et alj l2002f ) but simulations of disk 
galaxy satellite encounters indicate that interactions in 
which the orbit of the LMC intersects the disk during a 
close encounter i n the past may excite a strong tidally 
generated spiral ([Weinberg fc B litz 2006|). 

A close encounter with the LMC that brings it to 
half its current distance is expected to correspond to a 
quadrupolar tidal field that is about 8 times stronger 
since c ~ r~ 3 . Thus, the ratio qudai = c 2 ,mw/c2,lmC is 
a good indicator of the active strength of the halo tidal 
field on Galactic dynamical evolution. If qudai ~ 1 then 
the tidal field is comparable to the effect of the LMC 
on the Galaxy and so is quite weak while if qudai ~ 8, 
we might expect significant tidal torques on the Galaxy 
that may cause interesting evolution in the form of spi- 
rals, bars and warps. 




Fig. 1. — The distribution of the relative alignment between the 
inner and outer halo from the sample of cosmological dark halos. 
We measured the normalized moments of inertia for particles in 
the inner halo with r < r s and the outer halo for particles in the 
range r s < r < 2r s . Diagonalization of the inertia tensor deter- 
mines the direction of the principle axes. We calculate the angle 
between the minor axes of the inner and outer halos 59 and plot 
the distribution. As expected the inner and outer halos are closely 
aligned though the median of distribution around 20° implies a 
significant misalignment between in the inner and outer halos for 
at least half of the dark halos. 

Figure shows the distribution of values of the pa- 
rameter qudai measured for the sample of halos using 
two different outer radii, ro = r s and ro = 2r s , plot- 
ted against the axis ratio c/a. The halo axis ratios are 
determined using the normalized inertia tensor of halo 
particles within r < r s using an iterative method that 
gives a measure of the best fit e llipsoid of the distribu- 
tion (Dubinski & Carlbcrg 1991). There appears to be 
no strong correlation between the value of c 2 and halo 
shape since the scatter is quite large. The mean value 
of qudai is about 5 for r = 2r s compared to about 25 
for a cutoff at ro = r s . The implication then is that a 
Milky Way disk that is misaligned with the outer halo 
will experience a tidal field that is perhaps 5 to 25 times 
the strength of the tidal field of an LMC-type satellite, 
depending on the radial location where the misalignment 
begins. For halos with a larger concentration c, both r s 
and M 2 oo are smaller for the MW model but the values 
of c 2 vary within a factor of a few. This result is some- 
what surprising for it implies that the external tidal field 
acting on a disk galaxy from a triaxial halo is comparable 
in strength to the tidal field of a large satellite having a 
close encounter with a galactic disk. 

Finally, the disk will only be torqued if there is a sig- 
nificant misalignment with the halo independent of the 
strength of qudai ■ One might imagine that more flattened 
halos with the largest qudai are more closely aligned and 
thus have a weaker potential for torquing. We therefore 
looked for a correlation between qudai and the relative 
alignment of the inner and outer halo 89. We see no 
such effect (Fig. [3]) and conclude that significant per- 
turbations can arise from external torquing due to this 
misalignment. 

We now go on to demonstrate some of these effects in 
controlled N-body simulations. 
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Fig. 2. — Relative values of the external quadrupolar tidal field 
of triaxial dark halos measured by the parameter, C2, for matter 
beyond radius r s (red) and 2r s (blue). To give the value of C2 a 
physical meaning in the context of the Milky Way, we calculate C2 
for a satellite similar to the Large Magellanic Cloud with a mass, 
M = 10 10 Mq at a distance r = 50 kpc. We then calculate C2 
for each dark halo and rescale the model to a Milky Way galaxy 
model with mass M200 = 2 X 10 12 Mq and r s = 26 kpc. The ra- 
tio C2,mw I C2,LMC i s therefore a measure of the relative strength 
of the tidal field of the external quadrupole to a LMC satellite. 
It is known that the strength of the LMC tidal field is weak and 
has minimal dynamical effects but a field that is 5 to 10 times 
stronger can have a significant effect. The strength of the exter- 
nal quadrupole field is 5 to 25 times the strength of the LMC for 
matter beyond 1 to 2 scale radii suggesting interesting dynamical 
phenomenology of the disk as the result of triaxial halos. 




Fig. 3. — A scatter plot of the normalized strength of the exter- 
nal tidal field qudal an d the misalignment angle between the inner 
and outer halo 58 for 2200 dark halos. There is no correlation, sug- 
gesting that a significant fraction of galactic disks may experience 
relatively strong tidal torques from halo misalignment during their 
evolution. 



3. SIMULATIONS 

We perform a series of TV-body simulation s using a 
modi fied version of a parallelized treecode (|Dubinskil 
1996). The code has been changed to include an exter- 
nal quadrupolar potential field that is tumbling at a fixed 
pattern speed Q p . Forces on particles are the sum of the 
A^-body force and the external potential field. If parti- 
cles stray beyond a radius rn then the field is shut off. 
For the simulations described here, the radius is ro = 60 
kpc. A model of a galactic disk plus surrounding dark 
halo is simulated within this tidal field. We also make 
sure the model remains centered on the quadrupole field. 
This prevents the outer dark matter halo in the N-body 
model from being significantly distorted. 

The model galaxies in question represent the 
Milky Way and Androm eda (M31), as given by 
IWidrow fc Dubinskil (|2005[ ) (WD herein). These models 
both have exponential disks with Hernquist-type bulges 
embedded within a cuspy dark halo similar to an NFW 
profile. We use the most stable versions of these mod- 
els called MWb and M31a in WD but we refer to them 
simply as MW and M31 in this paper. The M31 model 
has been tested at a variety of resolutions and is known 
to b e stable against bar f ormation after run times of 10 
Gyr (jGauthier et al.ll2006| ). The MW model on the other 
hand is mildly unstable and develops a bar later on af- 
ter t = 6 Gyr as we shall see below. The particles are 
distributed with the following numbers in both the M 
31 and MW runs: AT disfe =l,500,000, N bu i ge =500,000 and 
A/j a / o =2,000,000. The model parameters are provided in 
Table 1. 

We also populate the disks with particles out to large 
radii to better investigate edge effects. The MW model 
extends to 10 radial scale lengths while the M31 model 
goes out to 7 scale lengths. We note that the scale radii of 
the MW and M31 halos are r s k 10 kpc and somewhat 
smaller than our putative MW halo in the discussion 
above. The parameters are fitted values for observed 
rotation curves and inferred brightness profiles of the two 
galaxies. The real inner halos of M31 and MW are likely 
contracted compared to the pure dark matter models. 
The mass model fits to the real data have scale radii 
that are smaller by a factor of two in accordance with 
expectations of the contraction of halos by dissipative 
processes during galaxy formation. 

We have found it convenient to select a set of simula- 
tion units that set the scale lengths and rotation veloc- 
ities to values near unity. We scale the simulation units 
to the physical units in the following way: 

• Gravitational constant G=l, 

• 1 simulation length unit = 4 kpc 

• 1 simulation velocity unit = 220 km s" 1 , 



• 1 simulation mass unit = 4.5 xlO 10 Mp 



0i 



• 1 simulation time unit = 17.7 Myr 
We express the general external tidal field as: 



m—l 



$ext{r,0,<j)) = ^2 ^2 Pim{cos 9) (aim cos m(f)+bi m sin 7 

l m=0 



(10) 
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TABLE 1 

Parameters for the Milky Way (MW) and M31 models used in our simulations (taken from Widrow & Dubinski, 2005) 



Model M d R d z d r s v s v b r b Rq 5Ro Q 
MW 3.3 2.81 0.44 8.8 345 435 0.88 30 1 1.3 at 2.5R d 
M31 7.7 5.58 0.60 12.9 337 461 1.83 30 1 1.25 at 2R d 

Col. 1: Model for galaxy, Col. 2: Disk mass (10 10 M Q ), Col. 3: Disk scale length (kpc), Col. 4: Disk scale height (kpc), Col. 5: 
NFW-halo scale radius (kpc), Col. 6: NFW-halo scale velocity (km s — 1 ), defined in Section 2, Col. 7: characteristic bulge velocity 
(km s _1 ), Col. 8: bulge scale length, Col. 9: disk truncation radius (kpc), Col. 10: truncation width (kpc), Col. 11: Toomre 
Q-parameter. 



The coefficients a\ m and 6/ m describe the strength of the 
field and are derived from Equation [5] If we assume that 
the quadrupole is aligned with the principle axes and 
is due to a plane-symmetric distribution then the bi m 
terms are zero. For these simulations we are only inter- 
ested in the I — 2 quadrupole terms. The pure external 
quadrupole field then becomes: 

<b ext {r,6,4>) = ^2 a2 m P2 m (cos6) cos(rri(t))r 2 . (11) 

We estimate a suitable tidal field using the following 
procedure. We first take an NFW profile that has been 
flattened into a perfect ellipsoidal shape with axis ratios 
qi — b/a and qi — c/a. We then determine the spheri- 
cally averaged density profile of this distribution and fit 
a spherical NFW profile to determine the scale radius, 
r s and mass M s as is done in cosmological simulations. 
We are interested in the tidal field due to material be- 
yond a radius, ro = 2r s , so we compute the coefficients 
of the quadrupolar tidal field for the material beyond 
this radius. As an example, we find that for a flattened 
NFW halo with qi = 0.6 and g 2 = 0.5 with a scale length 
r s = 26 kpc and M 2 oo = 2 x 10 12 M Q that the quadrupo- 
lar tidal field coefficients in simulation units are: 

• a 20 = 0.00 50, 

• a 2 i = 0.0, 

• a 22 = -0.000 1 5. 

• The quadrupole amplitude is c 2 = 0.00051. 

(If we used a more concentrated halo with c = 15 then 
r s = 16 kpc and M 20 o = 1-5 x 10 12 M leading to a field 
strength that would be about 3 times larger for r = 2r s .) 
The N-body galaxy models have a nearly spherical halo 
at large radii so we model the effect of the triaxial halo 
by this quadrupolar field. We assume that the inner halo 
has aligned with the disk and is essentially axisymmetric 
and remains aligned for the course of the simulation. As 
a comparison, the quadrupole coefficients for a LMC- 
like satellite of M — 10 10 M Q , at a distance of 50 kpc 
placed on the £-axis, is c 2 = 8.12 x 10~ 5 in simulation 
units so that qudai = 6.3 close to expectation value for 
the cosmological distribution of halos for material with 
r > 2r s .. 

To introduce a misalignment, the disk is tilted by 30° 
with respect to the field and the field is made to ro- 
tate about the z-axis with different pattern speeds. This 
choice of the tilt angle is arbitrary but is approximately 
the expected angle for at least 30% of dark halos, ac- 
cording to the distribution shown in Figure [TJ The 4 
pattern speeds that we choose are fl p = 0.11,0.22,0.44, 



and 0.88 km s -1 kpc" 1 corresponding to tumbling pe- 
riods of T = 56, 28, 14, and 7 Gyr respectively. We note 
that these are very slow compared to the pattern speeds 
of barred galaxies (fi p w 20 km s -1 kpc -1 ) and so the 
corotation radii are mainly larger than r 2 oo (The cir- 
cular frequency at r 2 oo in the NFW model is f2 2 oo = 
f20o/f20o = #o/100 w 0.7 km s -1 kpc -1 .) These pe- 
riods are approximately 4.0, 2.0, 1. and 0.5 times the 
age of the universe. We recall that iBailin &: Steinmeti 
(2004) find a log-normal distribution peaking at fl p = 
0.15 km s _1 kpc -1 , corresponding to a tumbling period 
of 44 Gyr. Thus, our simulations probe the fast tumbling 
tail of this distribution - the half of the distribution that 
is most likely to yield interesting dynamical effects. The 
quadrupole coefficients and pattern speeds for the MW 
and M31 runs are summarized in Table 2. 

We simulate each run for 8000 equal timesteps with 
St = 0.05 units, corresponding to about 7.1 Gyr in phys- 
ical units. We soften gravity with a Plummer softening 
kernel with radius of 40 pc for the stars and 100 pc for 
the dark matter. Energy is conserved typically to within 
0.5% and total angular momentum to within 1%. 

4. RESULTS 

Before exploring the results for externally torqued 
disks, we first examine the evolution of control models. 
Figure |4] shows the final state for the MW and M 31 
models in the absence of any external torquing potential. 
The disks remain co-planar and there are minimal signs 
of disk thickening and warping resulting from the ampli- 
fication of the Poisson noise in the N-body disk. There 
are some remnants of spiral structure in the disks which 
arise from swing amplification of the noise in the disk. 
This transient spiral structure heats the disk azimuthally 
and slowly fades away as the simulations progress. 

When we turn on the external quadrupole potential 
with the expected amplitudes from cosmological dark ha- 
los, we clearly see the various effects predicted for disk 
torquing. Figures [5] and O show the evolution of repre- 
sentative models of the MW and M 31 runs from 3 per- 
pendicular views. Associated videos 3 also show how the 
external torque causes the disk to precess and nutate like 
a tilted gyroscope on a table top. The evolution of the 
direction of the spin vector of the disks described by the 
normalized x and y components of the angular momen- 
tum quantifies this behavior (Fig. [7]) and the nutational 
frequency depends on the pattern speed of the forcing 
quadrupole potential. All runs indicate that the model 
disk precesses through an angle of 30-60° over 7 Gyr as it 
responds to the external quadrupolar field implying pre- 
cession periods of more than 40 Gyr or precession rates 

3 Animations are available at the website 
http: / / www.cita.utoronto.ca/ ~dubinski / warpmovies 
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TABLE 2 

Simulation parameters 



Model 


i2o(icr 4 ) 


a22(10~ 4 


Op 


T 


MW-con 


0.0 


0.0 


0.0 


oo 


MW-0 


5.0 


-1.5 


0.11 


56 


MW-1 


5.0 


-1.5 


0.22 


28 


MW-2 


5.0 


-1.5 


0.44 


11 


MW-3 


5.0 


-1.5 


0.88 


7 


MW-2a 


2.5 


-0.75 


0.44 


28 


MW-2b 


1.3 


-0.38 


0.44 


28 


MW-2c 


0.63 


-0.19 


0.44 


28 


M31-con 


0.0 


0.0 


0.0 


oo 


M31-0 


5.0 


-1.5 


0.11 


56 


M31-1 


5.0 


-1.5 


0.22 


28 


M31-2 


5.0 


-1.5 


0.44 


14 


M31-3 


5.0 


-1.5 


0.88 


7 



Col 1: model name, Col 2: component 020 of the external quadrupolar field (10 4 
quadrupolar field (10 — 4 Mqs - 2 ), Col 4: halo pattern speed Q p (km s -1 kpc — 



Mqs ), Col 3: component (122 of the external 
Col 5: tumbling period (Gyrs) 




Fig. 4. — The state of the Milky Way and M31 disks (left and right respectively) at the end of control simulations with no external tidal 
potential. While a very weak oval distortion can be noticed in the Milky Way disk at this time (f« 6.4 Gyrs), the M 31 disk does not show 
any signature of a bar. The warping of the disks is also absent and the disks remain in their initial plane. 



less than 0.15 rad Gyr -1 . We see below that these pre- 
cession rates are in accord with analytical estimates from 
a rigid disk model forced by an external potential. These 
precession rates are about an order of magnitude smaller 
than those found in mod els with a disk misaligned within 
a triaxial potential (e.g. JKui iken 199lt Ueon et al.ll2009h 
used in models of warping behavior that postulate that 
the main source of torque on a disk is due to he potential 
of the inner part of the triaxial halo. The phenomenol- 
ogy that we are exploring is therefore different and the 
resulting effects are slower and more subtle. So while 
the precession periods are long compared to the Hubble 
time, a disk can slew through and a significant angle over 
its lifetime and this dynamical evolution leads to a more 
gradual transient warping at the disk edge in a mecha- 
nism similar to those invoking the a ccretion of angular 
momentum (|Ostriker fc Binnevl ll989). We conclude that 
this must be an important process in disk dynamical evo- 
lution over the age of the universe. For real galaxies, we 
expect a range of halo tumbling pattern speeds, tidal 
field strengths and initial tilt angles so there should be 



also be a range of warping behavior in the population at 
large. 

4.1. Warps 

While the the disk stars within approximately 5 radial 
scale-lengths stay within the plane and act like a rigid 
body because of their self-gravity, stars near the edge of 
the disk begin to precess at different frequencies and be- 
gin to warp away from the disk. Our simulations indicate 
that the model disks develop strong warping early on in 
the runs and that such warping is sustained through the 
length of the run (about 7 Gyrs) (Fig. [8]). 

The disk demonstrates strong flocculent spiral pattern, 
to even the disk-edge where self-gravity is considered too 
weak to excite such structure. This is advanced as the 
handiwork of the external quadrupole since we do not no- 
tice prominent outer spiral patterns in the control runs. 
This is particularly noticeable in the MW runs that have 
a more exte nded disk. 

Following iLevine et alj (|2006t ). we quantify the disk 
warps that develop in the simulations by analyzing the 
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Fig. 5. — Snapshots of the Milky Way disk, at three different times from the four different runs near the beginning, middle and end of 
the simulation. The three views of the disk are presented from the point of views that are initially face-on (x — y top-left) and the two 
perpendicular edge-on views (x — z bottom- left; y — z top-right). The disk precesses in all cases and the tilting induces a transient warp in 
the outer regions of the disk that persists for about 7 Gyr. The tidal forcing also excites spiral structure in the outer disk. 



vertical deviation of the disk from a plane at different 
radii through the function, 

w(</>) = ^w m e im * ( 12 ) 



where <p is the particle azimuth. 

The inner, nearly co-planar disk is considered to be in 
the x — y plane. The particles inside the j th ring - of ra- 
dius Rj and width SR - are assigned the height Zj from 
the midplane with defined as the corresponding par- 
ticle azimuth. We first rotate into the frame of the inner 
disk before the analysis. We determine the midplane of 
the disk by computing the moment of inertia of the disk 
particles within 3 scale lengths where the disk is nearly 
planar and then diagonalizing the inertia tensor to give 
the orientation of the disk. We determine the coefficients 
W m using a Fourier analysis of the particle distribution 



through, 

The amplitude of the deviation is given by \W m \ and the 
phase angle through <f> — tan _1 [/m(W m )/ Re(W m )]/m. 
We can then construct functions of W m versus R to ex- 
amine disk warping. The function Wq(R) corresponds 
to 77i = bowl-shaped deviations, W\{R) corresponds to 
the usual m = 1 integral-sign shaped warps while W% (R) 
are second order "scalloped" warps. 

Figure [9] shows the results of this analysis for the 4 
models of both the MW and M31 system. The dominant 
vertical deviation is the m = I integral-signed warp. For 
the MW models, the warp turns upwards at R w f 5 kpc 
and continues to R = 30 kpc reaching about 5 kpc above 
the plane. The M 31 models shows a qualitatively similar 
behavior though the warping begins at the larger radius 
of R w 20 kpc but only reaches about 3 kpc above the 
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Fig. 6. — Final snapshots of two representative M 31 runs shown at t = 7.1 Gyr. The left panel corresponds to the three views of the 
M 31 disk (with the same layout as Fig. [5j at the end of the run M31-3 with the fastest forcing pattern speed. The panel on the right 
shows the end of run M31-0 with the slowest pattern speed. The warping of the disk is more subtle in this case but is still apparent. 




Fig. 7.— The evolution of spin axis of the MW and M31 models - MW-0, M31-0 (red); MW-1, M31-1 (green); MW-2, MW-3 (blue); and 
MW-3, M31-3 (black). This time evolution is represented here by the normalized components of the disk angular momentum vector L x /L 
and Ly/L. The circle corresponds to a fixed inclination of 8 = 30°. The tumbling axis of the quadrupole potential is perpendicular to the 
plane of this plot. The disk responds by precessing in the counter-clockwise direction as expected for this quadrupole potential. Nutations 
are clearly seen as the disk responds to the triaxial halo. The nutational period is generally smaller than the precession frequency. Even 
after the 7 Gyr time of the simulation the disk only precesses through about 60 degrees. 



plane at a radius of R = 40 kpc. The M31 model is more 
extended and massive than the MW model while being 
forced by the same external tidal field and so shows a 
smaller degree of warping behavior. 

We also examined the alignment of the warps by mea- 
suring the line of nodes (LON) from the phase angle de- 
termined from the vertical Fourier analysi s. We present 
the tip-LON plots for our models following iBriggsl (|1990l ) 
(Fig. [T0|) . The warps begin with a relatively straight line 
of nodes when they begin to move out of the plane at the 
disk edge but generally show a leading spiral curving out 
to ~ QRd following the observations from Brigg's anal- 
ysis. At larger radii, the warping becomes less coherent 



resulting mainly by the onset of differential precession of 
the stellar disk orbits and the modulated forcing of the 
tumbling external quadrupole potential. While we are 
only modeling a pure N-body system one expects that 
the mechanism will generate coherent warping in the gas 
as well at least out to 6 scale lengths. 

4.2. Bars 

The tidal distortion due to a galaxy interaction in a 
fly- by can trigger the ba r instability under some condi- 
tions (e.g. Noguchi 19871 : IGerin et al.lll990h . The effect 
is due to a transient external quadrupole that disturbs 
the stellar orbits in the center of the galaxy. Under our 
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Fig. 8. — A prominent warp develops by the end of the simulation in model M31-3 with a quadrupole field pattern speed of Q p 
0.88 km s- 1 kpc" 1 . 



hypothesis of a tumbling misaligned halo, there is also a 
changing external quadrupole potential that can perturb 
the disk and potentially trigger the bar. The physical 
situation is somewhat different than a impulsive galac- 
tic fly-by however, in that rate of change of the orien- 
tation is slow and the amplitude of the tidal potential 
is constant. Nevertheless, the Milky Way disk is found 
to develop a fast bar at about 3 Gyrs. Once formed, 
the bar is sustained and its pattern speed declines from 
an initial value of 50 to 35 km s _1 kpc -1 by the end 
of the run in accord with typical expectations for mod- 
els like this fFigfTT]) (e.g . . iDebattista k. Sellwo"odl ll998t 
lO'Neill fc Dubinskill2003l : iDubinski et al.ll2OO90 . 

Figure [12] shows the evolution of the strength of the bar 
that is triggered in the Milky Way for different models 
by the torquing provided by the external halo. We find a 
well-developed bar, of nearly similar strengths, forming 
by about 3 Gyrs in all the runs, indicating that the rate 
of tumbling of the halo is immaterial to the bar strength. 
However, fi p of the quadrupole appears to have a weak 
effect on the time of onset of the bar instability - this is 
hastened with increasing tumbling period T, i.e. falling 
pattern speed of the external halo, though this is not 
marked; the fractional difference in bar onset times be- 
tween the runs done with the slowest and fastest halo in 
the Milky Way model (MW - and MW - 3, respec- 
tively) is about 16%. This small difference in the bar 
onset time is also noticed in the snapshots from these 
two runs, shown in Figure [5l 

Instead, we find that it is the strength of the 
quadrupole of the external halo that is the crucial fac- 
tor in controlling the bar instability. In fact, a control 
run performed with a null external quadrupole leads to 
a much later (« 5.8 Gyrs) onset of a very weak bar in 
the Milky Way disk (maximum bar strength is less than 
0.4 times the bar strength reached in the other runs); see 
Figure SI Thus, we assure ourselves that the inherent bar 
unstable nature of the Milky Way model is very weak, 
albeit non-zero. Additionally, the control runs for both 



models result in significantly weaker spirality, especially 
at the edges of the disk and no warp. 

The M 31 disk is found to be unaffected by the halo 
torquing, in regard to bar formation. This is apparent 
from the pictures of the M 31 disk at the end of the 
simulations (Figure [6|) ; the result is found independent 
of the halo pattern speed used in the run. 

5. ANALYTIC RIGID DISK MODEL 

It is interesting to gauge the effects of the external 
quadrupolar term in the potential of a tumbling triax- 
ial dark halo in an analytical model using a rigid, thin 
exponential disk for comparison to the phenomenology 
in the simulations. We emphasize again that the inner 
disk and halo are tightly coupled through dynamical fric- 
tion and so remain aligned. The only possible source of 
torquing on the disk is a misaligned and possibly tum- 
bling outer halo. The assumption of a rigid disk is a rea- 
sonable assumption since the simulations show that the 
disk maintains coherence because of self-gravity within 
R < 5Rd even while it tips in response to the external 
tidal torque. We use the classical Euler equations of mo- 
tions for a rigid spinning disk. The dynamic variables 
are the standard Euler angle (9, <fi, ip) where 9 is the 
inclination of the disk, <f> is the longitude of ascending 
node and if) the angle about the disk's symmetry axis. 
In the presence of a time dependent torque, the disk will 
precess and nutate and so 9 and 4> will evolve in time. 
In these calculations, we set up the disk with an initial 
inclination of 9 = 30° and set the halo quadrupole po- 
tential tumbling about the z-axis with 9 = for direct 
comparison to the results of the simulations.. 

5.1. Euler Equations of Motion 

After Goldste in et alJ (|2002l ) (see p. 210 and equation 
5.52), the Lagrangian of the disk is 

£j.(02 + 02 gin 2 0) + ^(0 C os 6> + tp) 2 - V{6, tf>,t), 

(14) 



C 
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Fig. 9. — Disk warping behavior quantified by the functions Wo(r) (dotted lines), Wi(r) (solid lines), and Ws(r) (dashed lines) for 
the four MW and M31 runs with halo tumbling periods of P = 56, 28, 14, and 7 Gyr (black, red, green and blue) respectively at the last 
snapshot at t = 7 Gyr. Disks generically develop vertical deviations over the radial range of 5, 6, and 7 R<i (vertical dashed lines) and 
are dominated by the m = 1 integral-signed warp with stars reaching 2-5 kpc above the plane depending on the model. The quantitative 
behavior of these warps are similar to that seen in real galaxies including the Milky- Way and M 31. 



where I3 is the moment of inertia about the symmetry 
axis l\ = I2 are the moments measured about an axis in 
the plane of the disk. The potential energy V is deter- 
mined from the interaction between a rigid thin exponen- 
tial disk and a tumbling external quadrupole potential. 
The equations of motion can then be written as: 



dV 



h9 -htf sin0cos0 + S<£sin0+ — = 0(15) 

o9 

dV 

h4> sin 2 9 + 2I 1( f>9 sin 9 cos 9 - S9 sin 9 + — = 0,(16) 

0(f> 

where S = l3((j)COs9 + ip) (the disk angular momentum) 
is a constant of motion and V is the potential energy 
between the disk and the tumbling halo. A tilted disk will 
naturally precess. For example, with an axisymmetric 
potential V(9), the expected precession frequency can be 
derived from the first Euler equation by assuming that 9 



is constant. For small 



we expect 

1 dV 



fprec 



S sin 9 89 ^ 
(cf. iGoldstein et al. 2002) . If the potential also depends 
on <p and is rotating about the z-axis, we can expect a 
more complex evolution combining both precession and 
nutation. 

5.1.1. Direct Solution of the Euler Equations for a Rigid 
Exponential Disk 

We first consider the interaction of a rigid exponential 
disk with an external quadrupolar tidal field for com- 
parison to the results of the N-body simulations. We 
point out that this process is quite different from the 
disk precession examined in work that ha s treated a tri- 
axial halo as a rig id background potential ( Kuijkcn| |l991l : 
iJeon et ail I2009D . In those treatments, the precession 
rate is an order of magnitude larger due to the much 
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Fig. 10. — The inclination (tip) versus longitude or tip-LON plots for the warps in the MW disk models at the last snapshot of the 
simulations at t = 7 Gyr. The four models with tumbling periods of P = 56, 28, 14, and 7 Gyr are give by the black, red, green, and 
blue lines respectively. The dashed circles correspond to 3 and 6 degrees of tip respectively while the longitude is measured around the 
circles. The line of nodes of the warp is approximately straight when it begins to move out of the plane but becomes more erratic at 
higher inclinations. The line of nodes tends to curve in a leading spiral form similar to what is observed in real warps. The warps are less 
pronounced in the M31 models but still behave in a similar manner. 




t (Gyr) 



Fig. 11. — Pattern speed evolution for model MW-2 from t = 
4 Gyr. The disk precesses initially in response to the external 
torques from the quadrupolar field and at late times develops a 
bar instability. The bar pattern speed decays in the usual manner 
in response to resonant interactions with the dark halo. 

larger torques felt by the inner halo. However, the large 
torque in these simulations is artificially high since we ex- 
pect dynamical friction to align th e disk and inner halo 
withi n a few disk dynamical times (jDubinski fc Kuijkenl 
1995). In this study, we assume that only the outer halo 
can exert any torque and subsequently the strength of 
the torque is considerably smaller, typically by an order 
of magnitude. 

In Appendix A, we compute the time dependent in- 
teraction potential between a rigid exponential disk and 
an external quadrupole potential to plug into the Eulcr 
equations. The potential is: 

V(9, ct>; t) = MI d R 2 d { - (3 cos 2 6-1) 



Fig. 12. — The evolution of the bar strength A2 is shown here 
for the 4 Milky Way runs as well as for the control run done with 
the Milky way model. Here, the bar strength is measured in terms 
of the quadrupolar strength parameter of the bar that forms in 
the disk. The broken black line corresponds to the fastest halo 
{MW — 3), while the green curve refers to MW — 2, followed by 
the blue and red curves, which depict A2 evolution for MW — 1 
and MW — 0. The control run is represented by the solid black 
line. 

+3a 2 2cos[2(ry-0)]sin 2 6»}. (18) 

where Md and Rd are the mass and exponential scale ra- 
dius of the disk, 020 and 022 are the quadrupole terms for 
the external terms of a triaxial potential in equation 1111 
fl p is the tumbling frequency of the quadrupole potential 
assumed to be about the z-axis with 6 = and we as- 
sume that <j) — at t — 0. We can compute the moments 
of inertia I\ and I3 and the spin S directly from the N- 
body model and so fully specify the Euler equations for 
our case. 
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For the M31 model, we solve the Euler equations of 
motion numerically using the same initial state as the 
N-body model along with four different pattern speeds 
for the quadrupole potential. We use the values of 
I\ = SMdRj and I3 = 6MdR d expected for an expo- 
nential disk as well as the spin S (disk angular momen- 
tum) measured from the N-body models. The specific 
parameters for the M31 model in dimcnsionless units are 
M = 1.75, R d = 1.39 so that h = 10.17, h = 20.34 
and S = 5.14. Using equation [T71 we can estimate 
the precession period assuming that 022 time-averages 
to zero for sufficiently large tt p . We find a value of 
4>prec — —0.14 km s _1 kpc" 1 which is comparable in 
magnitude to our chosen pattern speeds but opposite in 
sign. This halo corresponds to our choice of qudai of 
about 6.3. We then integrate the Euler equations for ap- 
proximately 53 Gyr corresponding to a little more than a 
typical precession period to get a complete picture of the 
evolution of the disk for comparison to the N-body sim- 
ulations. Figure fT3l shows the evolution of the spin axis 
of the disk for comparison to the similar plots for the N- 
body simulation (Fig.[7|)- As in Fig. [3 Figure[T3lprcscnts 
this evolution by tracing the normalized, x component of 
the disk angular momentum against the y component of 
the same. 

The agreement in behavior for the four pattern speeds 
is excellent suggesting that the N-body disk is behaving 
in a similar manner to a rigid body over this time pe- 
riod. The nutation induced by the triaxiality is readily 
apparent again in this plot. We notice that the nutation 
period is shorter for a higher pattern speed. Thus, for 
larger pattern speeds, as the nutation period is small, 
we might expect a spiral galaxy to wobble several times 
over its lifetime if it is embedded in a rapidly tumbling 
dark halo. For the slower pattern speeds, the nutation 
period becomes comparable to the precession period and 
the resulting effect is simply a steady change of the disk 
orientation. Stars and gas on the outer edge of the disk 
that have weaker self-gravity are prone to precess differ- 
entially and this may be an additional origin of galactic 
disk warps. 

5.2. Resonant Interactions 

For the flattened halo discussed here the precession 
rate of disk is negative while the tumbling rate is posi- 
tive and in the same sense as the disk. While this case 
is the most reasonable, it is possible that the sign of 
the precession rate and halo tumbling rate might be the 
same and so resonances can occur when Q p = <j) pre c- In 
the firs t investigations of orb its in a rotating triaxial po- 
tential, iBinneyl (| 1978. 1981) discovered that under cer- 
tain conditions retrograde orbits in the rotating frame 
can be unstable leading to large ^-motions due to reso- 
nant interactions with the potential and suggested that 
this might b e a me chanism to produce galactic warps. 
iHeisler et al.l |l982) demonstrated the existence of Bin- 
ney's resona nce in orbital studies us ing an analytic triax- 
ial potential. iTremaine &: Yul (|2000f ) have even suggested 
that this resonance may provide a mechanism for creat- 
ing polar ring galaxies. Accreted gas and forming stars 
counter-rotating with respect to main disk may be levi- 
tated out of the plane if there is an underlying rotating 
triaxial potential. 
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Fig. 13. — The trace of the spin axis through the normalized 
components of the disk angular momentum vector L x /L and L y /L 
for the rigid body model of M31 forced by an external quadrupole 
potential. Colored lines correspond to the same pattern speeds 
applied to N-body models with M31-0 (red); M31-1 (green); M31- 
2 (blue); and M31-3 (black). The agreement with Fig. [7] is very 
good suggesting that self-gravitating disks behave like rigid bodies 
in this context. However, different pattern speeds lead to different 
nutational frequencies and so the edges of disks may flop around 
at different rates. 

So far the study in this paper has elucidated how self- 
consistent N-body disks evolve in tumbling triaxial po- 
tential rotating in the same direction as the disk. In the 
other case of a counter-rotating disk, it is interesting to 
see if Binney's resonance manifests itself in disk evolution 
in our models. 

Let us first consider the simple rigid disk model. Con- 
sider the Euler equations again for the case where 4> = 
and 9 — corresponding to steady precession at a 
fixed angle 9 at a constant rate. For this to be the 
case in general, we need to find the condition such that 
dV/d(j) = at all times. Differentiating equation [18] 
with respect to 4> an d setting it to zero implies that 
sin[2(O p t — (f))] =0. If 4> = 4> P rect then this implies a 
resonance when Q p = <f> pr ec- With this resonant condi- 
tion, we can use the first Euler equation to find <fi from 
the quadratic equation 

dV 

- hef 2 sin 9 cos 9 + Scj> sin 9 + — = (19) 

o9 

where is known from eauation lA25l (see Appendix 1). 
For a flat rotation curve with speed v c the disk angular 
momentum is given by S — 2MdRdV c - If we assume that 
4> prec is small then to a good approximation it is given 

by 

9 Rd 

4> P rec ~ -7: (020+2022)0036*. (20) 

2 V c 

For our model of M31, with our estimate of 020 and 022 
for the coefficients of an external quadrupole potential 
and an initial disk inclination of 9 = 30°, the resonance 
occurs when n p r . es = 4> prec = —0.055 km s _1 kpc -1 . 
Thus we see that the size of this typical resonant preces- 
sion frequency is comparable to the expected tumbling 
frequencies of dark halos. 
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Fig. 14. — The trace of the spin axis measured by the normalized 
components of the disk angular momentum vector for the rigid 
body model of M31 forced by an external quadrupole potential near 
the resonant pattern speed. Colored lines correspond to the pattern 
speeds U p = U p , re s (red); U p = Q p , res /2 (green); U p = 2f2 f> , res 
(blue); and Q p = 4£7p ire3 (black). The resonance describes a region 
of unstable behavior where the disk can undergo large changes in 
inclination with irregular precession and nutation. 

We can examine the behavior of the disk preces- 
sional evolution near this resonance. In Fig. [TJ] we 
represent the disk spin axis evolution to bring out 
the unstable behavior of rigid models with £l p set to 
^p,res/2, ^p 5 res ; ^^p,res and 4i^p .res- The radically dif- 
ferent near-resonance disk behavior noted in this figure, 
marks the position of a region of instability. Here the disk 
inclination can change by large amount and the preces- 
sion and nutational behavior are quite irregular. 

We also consider an additional M31 N-body model 
for comparison to the rigid disk model. We evolve it 
in the same way as previous models but here examine 
a counter-rotating external quadrupole potential with 
Q p = —0.11 km s _1 kpc -1 to see if the irregular be- 
havior that is seen in the simple rigid disk model is re- 
produced in a self-consistent N-body model. Figure [15] 
shows a comparison of the time evolution of the polar 
angles (9, </>) representing the direction of the disk angu- 
lar momentum vector for the rigid disk and N-body disk 
models. (Here (f> is offset by 90 degrees from the usual 
definition as the angle of the line of nodes in the Eulcr 
equation.) For the rigid disk, we adjust the initial incli- 
nation slightly to 9 = 34° and find excellent agreement 
with the N-body disk evolution. We have integrated for 
the unusually long time of 25 Gyr to see the full tilt evo- 
lution. However, the amplitude of the quadrupole poten- 
tial could easily be twice as large based on the observed 
variance of qudai in cosmological halos and so this behav- 
ior would occur in half the time since from equation [201 
we see that the precession rate is directly proportional to 
the external quadrupole amplitude. It is remarkable that 
in this counter-rotating case, the disk inclination grows 
by 90° before reversing course. 

While the Binney resonance is derived in the context 
of closed collisionless orbits in rotating triaxial poten- 
tials, a version of it also appears to be operating col- 
lectively in a self-gravitating disk in the same context. 



Fig. 15. — The time evolution of the polar angles (8, 0) of the 
direction of the disk angular momentum vector for an M 31 rigid 
disk model (smooth curves) and N-body model (points). In both 
cases, the external quadrupole potential is counter-rotating with 
Q p = —0.11 km s — 1 kpc -1 , a value that is twice the resonant 
frequency for this model. The agreement between the rigid disk 
model and N-body simulation is excellent when we adjust the initial 
inclination to 8 = 34° for the rigid disk case. In this case, the disk 
inclination increases by almost 90° before reversing direction. 

Our idealized model of a N-body disk forced by a rotat- 
ing external quadrupole potential tips through 90° when 
counter-rotating tumbling frequencies are near the reso- 
nant value. The case of counter-rotation is probably rare 
but not impossible given the existe nce of some disks w ith 
counter-rotating components (e.g.. lRubin et al.lll992t ). It 
may also help e xplain the phen o meno n of polar rings 
as suggested by iTremaine fc Yul (|2000h . Accreted gas 
that is counter-rotating with respect to the tumbling halo 
might rise out of the halo principal plane within a Hub- 
ble time given a large enough quadrupole amplitude and 
suitable tumbling frequency. One caveat to be aware of 
is that given the long timescales of this process, dynam- 
ical friction may begin to align even the outer parts of 
the halo so the resulting torque may diminish with time 
and so weaken the effect. 

6. CONCLUSIONS 

In this paper we discuss the effect of the tidal field of 
a tumbling and triaxial external halo on the dynamical 
evolution of a galactic disk. We have shown that the 
disk precesses in response to its misalignment with the 
external halo and the gradual re-alignment of the disk 
with the external tumbling halo causes it to warp. The 
disk is also seen to develop strong spirality to its edge 
in the potential of the external halo that it sits in. It is 
also suggested that if a disk feels a significant external 
quadrupole that is essentially static or at most slowly 
tumbling, it might trigger a bar instability, late in the 
disk. 

The conclusions achieved in this work are based on 
analytical calculations that are backed up by N-body 
simulations of the Milky Way and M 31 models, sub- 
ject to the quadrupolar tug of an external tumbling and 
triaxial halo. The quadrupolar strength is judged via 
cosmological simulations, while a range of halo tumbling 
periods are scanned through in the simulations. We ad- 
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vance this dynamical mechanism as an important cause 
of non-axisymmetric structures to set in galactic disks. 
An important result of our work is that the prodding of 
the disk by a external quadrupole may not always pro- 
duce bars in the disks. Thus, the M 31 disk was found to 
be robust to the bar instability, though the Milky Way 
disk was found to be much more susceptible to the bar 
being triggered by this mechanism. However, both disks 
were found to warp. The dispersion in cosmological halo 
properties imply that the external quadrupole is rather 
weak in some systems and so may not be effective in 
driving disk evolution. 

Our work elucidates the importance of the quadrupolar 
term in the potential of the external halo, particularly, in 



terms of the trigerring of the bar and warp instabilities. 
Given the viability of such tidal prodding of the disk by 
the external halo, detailed investigation of this mecha- 
nism in galaxy formation simulations is suggested. 
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APPENDIX 



We calculate the interaction potential between an exponential disk and a time-dependent external quadrupole poten- 
tial rotating with pattern speed Q p . If we first assume the potential is static and aligned with the body axes (x' , y' , z') 
then we can write the potential from equation II II in Cartesian coordinates as: 

$ e * t (z',y', z) = Wx' 2 + b 2 y' 2 + b 3 z' 2 (Al) 

with 

x' = r cos (j> sin 8 (A2) 

y' = r sin <f> sin 8 (A3) 

z' = r cos 9 (A4) 

where one can show easily that the coefficients b±, b 2 and 63 are defined in terms of the harmonic expansion in 
equation QT] through: 

bi = -a 20 /2 + 3a 22 (A5) 

62 = -020/2 - 3a 2 2 (A6) 

63 = a 20 (A7) 

If the potential is rotating counterclockwise about the z'-axis with pattern speed f2 p then we can introduce the time- 
dependence of the potential into the coordinates (x' , y , z') using inertial frame coordinates (x, y, z) where 

x —x cos Q p t — y sin Q p t (A8) 

y' = x sin Q p t + y cos Q, p t (A9) 

z' = z (A10) 



so that potential becomes 
with the coefficients 



<S> ext (x,y,z;t) = dxififx 2 + d 2 (t)y 2 + d 3 {t)xy + d iZ 2 (All) 



di(t) = -a 20 /2 + 3a 22 cos2Q p t (A12) 

d 2 (t) — —a 2 o/2 — 3a 22 cos 2n p t (A13) 

d 3 (t) =6a 22 sin2n p t (A14) 

6^4 =020 (A15) 

We can thus derive the time-dependent interaction potential between a tilted exponential disk and this tumbling 
quadrupole potential. The parametric equation of a ring of radius R with inclination angle 6 and line of nodes <f> is 
given by: 

x = i?(cos0cos ip — sin cos sin -0) (A16) 

y = i?(sin cj) cos tp + cos (f> cos 8 sin ip) (A17) 

z = R sin 8 sin ip (A18) 

where tp is the angle measured along the ring. The mass element of the ring is given by 

dM = Z(R)R dR dip (A19) 

so we can derive the interaction potential by integrating over tp and R through 

V(8,(f>;t) = J d^j J Z(R)RdR$ ext {x,y,z) (A20) 
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If we assume an exponential disk with 



W) = ^e-W (A21) 



and substitute the expressions for (x, y, z) in equations IA18I and do the integrals we arrive at the intermediate form of 
the interaction potential: 

V{9,(f>;t) = 3M d B%[di{t)(l - sin 2 sin 2 6>) + d 2 (t)(l - cos 2 sin 2 6>) + +d 3 (t) (sin 2</> sin 2 60 /2 + d 4 sin 2 6>] (A22) 
This can be reduced to the following simpler form in terms of the original quadrupole coefficients 020 and 0,22'- 

V(9, 0; t) = 3M d R 2 d {-^(3 cos 2 9-1) + 3a 22 cos[2(fl p t - <j>)] sin 2 6»| (A23) 

The derivatives of the interaction potential used in the Euler equations are then: 

dV 

— = 9M d R 2 d sin 9 cos 9 {a 2 a + 2a 22 cos[2(fl p t - cp)} } (A24) 
06 

dV 

— = 18M d R 2 d a 22 sin[2(fl p i - (j>)} sin 2 9 (A25) 
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